      !=======================================================================
      !> @file extfunc.f90
      !> @brief User provided functions and subroutines module
      !> @authors Ary Rodriguez-Gonzalez & Alejandro Esquivel
      !> @date 23/Jun/2011
      !
      ! Copyright (c) 2011 Ary Rodriguez-Gonzalez & Alejandro Esquivel
      !
      ! This file is part of AGA-V1.
      !
      ! AGA-V1 is free software; you can redistribute it and/or modify
      ! it under the terms of the GNU General Public License as published by
      ! the Free Software Foundation; either version 3 of the License, or
      ! (at your option) any later version.
      !
      ! This program is distributed in the hope that it will be useful,
      ! but WITHOUT ANY WARRANTY; without even the implied warranty of
      ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
      ! GNU General Public License for more details.
      ! You should have received a copy of the GNU General Public License
      ! along with this program.  If not, see http://www.gnu.org/licenses/.
      !=======================================================================
      !> @brief Number of Gaussian components
      !> @details Short module to implement the number of Gaussian components
      !! used for fit (employed for HH34) as a variable of easy global access.
      module numgaus
      !> @param ngaus : number of Gaussian components to be fit
        integer :: ngaus
      end module numgaus
      !-----------------------------------------------------------------------
      !> @brief External functions module 
      !> @details Here the user needs to provide a
      !> target (tarfunc) and a merit function (meritfunc). We recommend to
      !! preserve the names and interface in order to avoid changing code 
      !! elsewhere.
      module exfunc
      contains
      !---------------------------------------------------------------
      !> @brief Target Function
      !> @details User provided target function. \n
      !! In case some observational or experimental data is being fit,
      !! this is the function adjusted.\n
      !! For the example implemented here, finding a maxima of a given 
      !! function, the evaluation of the level of fitness can be made directly
      !! inside the merit function, thus the target function  returns a zero
      !! (and it is actually not called). Three examples of target functions 
      !! can be found commented below: several Gaussian components, one
      !! or two Lorentzians, and a PNe luminosity function.
      !> @param[in] x : Independent variable (abscissa)
      !> @param[in] npar : Number of independent parameters
      !> @param[in] q(npar) : Vector of independent parameters
      !> @return Value of the target function
        real function tarfunc(x,npar,q)
          !
          use numgaus
          !
          implicit none
          integer, intent(in) :: npar
          real, dimension(npar), intent(in):: q
          real, intent(in) :: x
      !    real, sum0,dintg,intg
      !    integer, NumPNs
      !
      ! Funcion Gausianas
      !    select case (ngaus)
      !    case (1)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))
      !    case (2)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
      !              +q(4)*exp(-(x-q(5))**2/q(6))+q(7)
      !    case (3)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
      !              +q(4)*exp(-(x-q(5))**2/q(6))&
      !              +q(7)*exp(-(x-q(8))**2/q(9))
      !    case (4)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
      !             +q(4)*exp(-(x-q(5))**2/q(6))&
      !              +q(7)*exp(-(x-q(8))**2/q(9))&
      !              +q(10)*exp(-(x-q(11))**2/q(12))
      !    case (5)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
      !              +q(4)*exp(-(x-q(5))**2/q(6))&
      !              +q(7)*exp(-(x-q(8))**2/q(9))&
      !              +q(10)*exp(-(x-q(11))**2/q(12))&
      !              +q(13)*exp(-(x-q(14))**2/q(15))
      !    case (6)
      !       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
      !              +q(4)*exp(-(x-q(5))**2/q(6))&
      !              +q(7)*exp(-(x-q(8))**2/q(9))&
      !              +q(10)*exp(-(x-q(11))**2/q(12))&
      !              +q(13)*exp(-(x-q(14))**2/q(15))&
      !              +q(16)*exp(-(x-q(17))**2/q(18))
      !    end select
      
      ! Funcion Lorentzians

         ! real, parameter :: pi=3.141592
          ! ngaus=1
         ! select case (ngaus)
         !    
         ! case (1)
         !    tarfunc=q(1)/(pi*q(2)*(1.+((x+q(3)))/q(2)))
         ! case (2)
         !    tarfunc=q(1)/(pi*q(2)*(1.+((x+q(3)))/q(2)))&
         !           +q(4)/(pi*q(5)*(1.+((x+q(6)))/q(5)))
         ! end select
      !
      ! Planetary Nebulae Luminosity Function
      !    tarfunc=q(1)*exp(0.307*x)*(1.-exp(3.*(q(2)-x)))
      !
      !    NumPNs=100
      !    summ0=0.
      !    Do j=0,NumPNs-1 
      !       summ0=3.*exp(3*q(1))/(exp(3.*q(1))-exp(3.*A(0,j)))+summ0
      !    EndDo
      !    dintg=1.114*(exp(3.*q(1)-2.693*ml)-exp(0.307*q(1)))    
      !
      !    tarfunc=0.

	!Prueba y= A*x + B

	tarfunc=q(1)*x+q(2)
          return
        end function tarfunc
      !-----------------------------------------------------------------------
      !> @brief Merit Function
      !> @details User provided merit function. \n
      !! This function computes the fitness level of each individual. \n
      !! The program assumes that a low value has the highest level of fitness
      !! (since its primary application is to fit functions with a merit
      !! function that is an error measure, e.g. @f$\chi^2@f$). If one needs
      !! to find a maximum of a function, one can use its reciprocal as merit
      !! function.\n
      !! The interface of the function is designed to adjust a target function
      !! to a set of observational/experimental data, there some arguments
      !! become dummy if the merit function does not need to be
      !! called.
      !> @param[in] xobs(ndata) : Vector of data abscissae
      !> @param[in] yobs(ndata) : Vector of data ordinates
      !> @param[in] ysigma(ndata) : Vector of uncertainties in the ordinate
      !> @param[in] q(npar) : Vector of parameters to be adjusted
      !> @param[in] npar : Number of parameters adjusted
      !> @param[in] ndata : Number of data points
      !> @return Value of the merit function (lowest value=fittest individual)
          real function meritfunc(xobs,yobs,ysigma,q,npar,ndata)
          implicit none
          real, dimension(ndata), intent(in) :: xobs,yobs,ysigma
          real, dimension(npar), intent(in) :: q
          real :: chi2,yfunc,pi,a,b,sig,x0,y0,r
          integer, intent(in):: npar,ndata
          integer :: i,n
      !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !                 Chi^2 fitting                                    !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          chi2=0.
          do i=1,ndata
             yfunc=tarfunc(xobs(i),npar,q)
             chi2=chi2+((yobs(i)-yfunc)/(ysigma(i)))**2
          enddo
          meritfunc=chi2
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !                       TEST FUNCTION 1                            !
      !       f(x,y)=[16x(1-x)y(1-y)sin(n*pi*x)sin(n*pi*y)]^2            !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !    pi=3.141592
      !    n=9
      !    meritfunc=16.*q(1)*(1.-q(1))*q(2)*(1.-q(2))*sin(pi*n*q(1))*sin(pi*n*q(2))
      !    meritfunc=1.-meritfunc
      !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !                      TEST FUNCTION 2                             !
      !                    r=sqrt((x-x0)^2+(y-y0)^2)                     !
      !                  f(r)=cos^a(b*pi*r)*exp(-r^2/(2*sig^2))           !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !    a=2.
      !    b=3.
      !    sig=2.
      !    x0=0.25
      !    y0=0.25
      !    pi=3.141592
      !    r=sqrt((q(1)-x0)**2+(q(2)-y0)**2)
      !    meritfunc=(cos(b*pi*r))**a*exp(-r**2/(2.*sig))
      !    meritfunc=1.-meritfunc
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
      !                      TEST FUNCTION 3                             !     
      !       f(x,y)=(1+cos(12sqrt(x^2+y^2)))/(0.5(x^2+y^2)+2)           ! 
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       
      !    meritfunc=(1+cos(12*sqrt(q(1)**2+q(2)**2)))/(0.5*(q(1)**2+q(2)**2)+2)
      !    meritfunc=1.-meritfunc
      !
      !
          return
        end function meritfunc
      !-----------------------------------------------------------------------
      end module exfunc
      !-----------------------------------------------------------------------
      !=======================================================================
      !> @brief Input from data and options file.
      !> @details contains routines to read the input data and run-time
      !! parameters.
      module input
      contains
      !> @brief Observational/experimental data input
      !> @details Reads experimental/observational data from file. Data are
      !! ordered by columns (see example in ObsData.dat)
      !> @param[out] xobs(ndata) : vector of abscissae
      !> @param[out] yobs(ndata) : vector of ordinates
      !> @param[out] ysigma(ndata) : vector of uncertainties of the ordinates
      !> @param[out] ndata : number of experimental points.
        subroutine ReadData(xobs,yobs,ysigma,ndata)
          implicit none
          real, dimension(:), allocatable, intent(out) :: xobs,yobs,ysigma
          integer, intent(out) :: ndata
          !
          real, dimension(:), allocatable :: x1,y1,ys1
          integer :: i, status
          integer, parameter :: n=10000000
          !
          !   read the input file (3 columns, x,y,sigma_y, and determine the
          !   number of lines it contains
          open(1,file='ObsData.dat',status='old',iostat=status)
          allocate(x1(n),y1(n),ys1(n))
          ndata=0
          do i=1,n
             read(1,*,iostat=status) x1(i),y1(i),ys1(i)
             if (status /= 0) exit
             ndata=ndata+1
          enddo
          close(1)
          !
          !   allocate the arrays (size=number of lines in input file)
          allocate(xobs(ndata),yobs(ndata),ysigma(ndata))
          !
          !   fill the arrays with the temporary buffers and free memory
          xobs  (:) = x1 (1:ndata)
          yobs  (:) = y1 (1:ndata)
          ysigma(:) = ys1(1:ndata)
          deallocate(x1,y1,ys1)
          !
          return
        end subroutine ReadData
       !---------------------------------------------------------------------
      !> @brief  Reads runtime parameters from file and initializes arrays
      !> @details  Reads from file (InitCond.dat) all the runtime input 
      !! parameters, initializes the values of the hyperboxes, and allocates
      !! memory for all the parents and population.
      !> @param[out] nsteps : Number of iterations before the hyperboxes are
      !!  restored to their initial size
      !> @param[out] nruns : Number of times that the hyperboxes are returned
      !! to their initial size
      !> @param[out] nparent : Number of parents
      !> @param[out] nind : Number of individuals
      !> @param[out] nerr : Number of realizations performed to obtain an
      !! estimate of the errors
      !> @param[out] deltaconst : if equal to 1 the reduction of the hyperboxes
      !! is at a constant rate (see reduc), otherwise hyperboxes are reduced 
      !! with the standard deviation of each parameter.
      !> @param[out] nmig : turns on/off migration. If its value is 0, 
      !! migration is disabled, otherwise is enabled
      !> @param[out] ngaussdev :  : if value equal to 1, offsprinf is obtained
      !! with a Gaussian distribution, otherwise they are obtained with
      !! uniform random deviates.
      !> @param[out] igaus :  ????? Quiza deberia pasarse al modulo de las 
      !! gaussianas
      !> @param[out] npar :  number of free parameters
      !> @param[out] gensize : size of the new generation from each parent.
      !! Since the parent is kept gensize=num of sons + 1
      !> @param[out] reduc : reducing factor, if a constant speed is chosen the
      !! hyper-box is reduced according to: 
      !> @param[out] beta : Tolerance value
      !> @param[out] qbounds(2,npar) : limits of the original (user provided)
      !! hyperbox
      !> @param[out] q0(npar) : Initial guess, computed as the position at the 
      !! center of the initial hyperbox
      !> @param[out] delta0(npar) : Initial window half-size
      !> @param[out] fit(nind) : Array to store the fitness level of the
      !! population
      !> @param[out] arrmain(npar,nind) : Array for the entire population
      !> @param[out] arrpar(npar,nparent) : Array for the parents only
      !> @param[out] arrgen(npar,gensize) : Array to store a single generation
      !! of each parent
        subroutine InitCond(nsteps,nruns,nparent,nind,nerr, deltaconst, nmig, &
                                                ngaussdev, igaus, npar, gensize, reduc, beta,     &
                                                qbounds,q0,delta0,fit,arrmain,arrpar,arrgen )
          implicit none
          integer, intent(out):: npar,nsteps,nruns,nparent,nind,gensize, &
                                  deltaconst, nmig,ngaussdev,nerr,igaus
          real, intent(out)   :: reduc,beta
          real, dimension(:),   allocatable, intent(out):: q0,delta0,fit
          real, dimension(:,:), allocatable, intent(out):: arrmain,arrpar,arrgen,qbounds
          real, dimension(:,:), allocatable :: aux
          !
              integer :: i,status
          integer, parameter :: n=10000
          !
          open(1,file='InitCond.dat',status='old',iostat=status)
          allocate(aux(2,1000))
          !
          !
      2   continue
          read(1,*,err=2) nsteps,nruns,nparent,nind,nerr
          !
      3   continue
          read(1,*,err=3) deltaconst, nmig ,ngaussdev, igaus
          !
      4   continue
          read(1,*,err=4) reduc,beta
          !
      5      continue
          read(1,*,err=5)aux(1,1),aux(2,1)
          !
          npar=1
          do i=2,n
             read(1,*,iostat=status)aux(1,i),aux(2,i)
             if (status/=0)exit
             npar=npar+1
          enddo
          !
          close(1)
          ! 
          allocate(q0(npar),delta0(npar))
          allocate(qbounds(2,npar))
          qbounds(1,:)=aux(1,1:npar)
          qbounds(2,:)=aux(2,1:npar)
          q0(:)=(qbounds(1,1:npar)+qbounds(2,1:npar))/2.
          delta0(:)=(qbounds(2,1:npar)-qbounds(1,1:npar))/2.
              !
              deallocate(aux)
          !
          gensize=nind/nparent
          allocate(fit(nind), arrmain(npar,nind),arrpar(npar,nparent),         &
                   arrgen(npar,gensize))
          !
          return
        end subroutine InitCond
        !--------------------------------------------------------------------
      end module input
      !----------------------------------------------------------------------
